# Seeing how emissions have changed between regions 2007-2015-2019
library(pacman)
p_load(
  here, fst, data.table, purrr, collapse, janitor, sf, tigris,
  dplyr, ggplot2, stringr, readxl
)


# Reading unit info dt with region names 
unit_info_dt = 
  read.fst(
    here("Data/electricity-generation/unit-info-dt.fst"),
    as.data.table = TRUE
  )[,.(orispl_code, unitid, nerc_adj)]

# Loading the electricity generation data
elec_gen_dt=
  read.fst(
    path = here("Data/electricity-generation/elec-gen-dt.fst"),
    as.data.table = TRUE
  )|> 
  setkey(orispl_code, unitid, utc_time) 
# Reading in pm2.5 emissions rates 
pm25_emissions_dt = fread(
  here('Data/electricity-generation/output/emissions-rate-avg.csv')
)
# Now for marginal damages
fips_dt = read_xlsx(
  path = here('Data/Marginal Damages (2011) from Holland, Mansur, Muller, Yates AER forthcoming.xlsx'),
  sheet = 'ground-level MD per ton'
) %>%  
  data.table() %>% 
  .[,.( # Miami-Dade changed their fips code
    fips = fcase(
      fips == '12025', '12086',
      fips != '12025', str_pad(fips, 5,'left','0')
    )
  )]
md_dt = 
  rbind(
    cbind(
      fips_dt,
      fread(here("Data/AP3/JointFolder/md_L_2014_CS.csv"))[,stack_height := 'low']
    ),
    cbind(
      fips_dt,
      fread(here("Data/AP3/JointFolder/md_M_2014_CS.csv"))[,stack_height := 'medium']
    )
  ) |> 
  setnames(new = c('fips','nh3','nox','pm25','so2','voc','stack_height')) 

# Calculating average damages by nerc region
avg_gen_dt = 
  elec_gen_dt[,.(
    tot_gload_mwh = sum(gload_mwh, na.rm = TRUE),
    tot_so2_tons = sum(so2_mass_lbs/2000, na.rm = TRUE),
    tot_nox_tons = sum(nox_mass_lbs/2000, na.rm = TRUE),
    tot_co2_tons = sum(co2_mass_tons, na.rm = TRUE)
  ), keyby = .(
    orispl_code, unitid, 
    nerc_adj, fips,
    stack_height = fcase(
      stack_height < 250, 'low',
      stack_height >= 250, 'medium'#, & stack_height <= 500,
      #stack_height > 500, 'tall'
    ))
  ] |>
  merge(
    pm25_emissions_dt[,.(orispl_code, unitid, tot_pm25_tons = pm25_tons_per_mw*tot_gload_mwh)],
    by = c('orispl_code','unitid'),
    all.x = TRUE
  ) |>
  merge(
    md_dt, 
    by = c('fips','stack_height'),
    all.x =TRUE
  )

avg_gen_dt[,.(
  avg_damage_per_kwh = sum(
    tot_so2_tons*so2 + tot_nox_tons*nox + tot_co2_tons*51*0.7993 + tot_pm25_tons*pm25,
    na.rm = TRUE
  )/(sum(tot_gload_mwh)*1000)
), keyby = nerc_adj][order(avg_damage_per_kwh)]



# Reading CEMS data for each year 
cems_dt =
  rbind(
    read_fst(
      here('Data/electricity-generation/cems/raw_cems_2007.fst'),
      as.data.table = TRUE
    )[,.(
      year = 2007,
      state, orispl_code, unitid, op_date, op_hour, op_time, 
      gload_mw = gload, sload_1000lb_hr = sload, 
      so2_mass_lbs = so2_mass, nox_mass_lbs = nox_mass, co2_mass_tons = co2_mass
    )],
    read_fst(
      here('Data/electricity-generation/cems/raw_cems_2019.fst'),
      as.data.table = TRUE
    )[,.(
      year = 2019,
      state, orispl_code, unitid, op_date, op_hour, op_time, 
      gload_mw, sload_1000lb_hr, 
      so2_mass_lbs, nox_mass_lbs, co2_mass_tons
    )] 
  )

# Checking state changes 
state_year_dt = 
  cems_dt[,.(
    tot_gload = sum(gload_mw*op_time, na.rm = TRUE),
    tot_so2 = sum(so2_mass_lbs, na.rm = TRUE),
    tot_nox = sum(nox_mass_lbs, na.rm = TRUE),
    tot_co2 = sum(co2_mass_tons, na.rm = TRUE)
    ), keyby = .(state, year)
  ][,':='(
    pct_tot_gload = tot_gload/sum(tot_gload),
    pct_tot_so2 = tot_so2/sum(tot_so2),
    pct_tot_nox = tot_nox/sum(tot_nox),
    pct_tot_co2 = tot_co2/sum(tot_co2)
  ), by = year]

state_diff_dt = 
  state_year_dt |>
  gby(state) |>
  fdiff() %>% 
  .[!is.na(year)] %>% 
  .[,year:=NULL] 

# Now getting shape files
states_sf = 
  states(cb = TRUE) |>
  filter( # Limiting to continental US
    !(STATEFP %in% c("02","15","60","66","69","72","78"))
  ) |>
  clean_names() |>
  st_transform(crs = 2163) |>
  left_join(
    state_diff_dt, 
    by =c('stusps'='state')
  ) 

theme_set(theme_minimal())

ggplot(data = states_sf, aes(fill = tot_so2)) +
  geom_sf(color = NA) +
  scale_fill_viridis_c('Change in SO2', direction = -1) 
ggplot(data = states_sf, aes(fill = tot_nox)) +
  geom_sf(color = NA) +
  scale_fill_viridis_c('Change in SO2', direction = -1) 










